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Abstract 



We show that a recently discovered fourth order symplectic algorithm, 
which requires one evaluation of force gradient in addition to three evalua- 
tions of the force, when iterated to higher order, yielded algorithms that are 
far superior to similarly iterated higher order algorithms based on the stan- 
dard Forest-Ruth algorithm. We gauge the accuracy of each algorithm by 
comparing the step-size independent error functions associated with energy 
conservation and the rotation of the Laplace-Runge-Lenz vector when solv- 
ing a highly eccentric Kepler problem. For orders 6, 8, 10 and 12, the new 
algorithms are approximately a factor of 10 3 , 10 4 , 10 4 and 10 5 better. 
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I. INTRODUCTION 



Symplectic algorithms JT|J^] for solving classical dynamical problems exactly conserve 
all Poincare invariants. For periodic orbits, the errors in energy conservation are bounded 
and periodic. This is in sharp contrast to Runge-Kutta type algorithms, whose energy error 
increases linearly with integration time, even for periodic orbits . Thus, symplectic algo- 
rithms are ideal for long time integration of equations of motion in problems of astrophysical 
interest M. For long time integrations, higher order algorithms are desirable because they 
permit the use of larger time steps. Symplectic algorithms are also advantageous in that 
higher order algorithms can be systematically generated from any low, even order algorithm 
In this work, we will show that higher order algorithms generated by a fourth order 
force gradient symplectic algorithm |J, have energy errors that are several orders of mag- 
nitude smaller than existing symplectic algorithms of the same order. For completeness, 
we will briefly summarize the operator derivation of symplectic algorithms and their higher 
order construction in Section II. While the materials in this section are not new, we be- 
lieve that we have restated Creutz' and Gockach's || triplet construction of higher order 
algorithms in its most transparent setting. In sections III and IV we recall force gradient 
algorithms and discuss two distinct ways of gauging the errors of an algorithm when solving 
the Kepler problem. We present our results and conclusions in Sections V and VI. 



II. OPERATOR FACTORIZATION AND HIGHER ORDER CONSTRUCTIONS 



After a tortuous start |10|,[11| , symplectic algorithms can be derived most simply on the 
basis of operator decomposition or factorization. For any dynamical variable W(qi,pi), its 
time evolution is given by the Poisson bracket, 

i,.„ , „„„, ^ f dWdH dWdH, 
If the Hamiltonian is of the form, 
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BM = \ y Erf + v{{*}), (2) 



the evolution equation (|1|) can be written as an operator equation 



dt ^ v dqi dp, 

with formal solution 



dW d ^ d \ , , 



W{t) = e*( T+y ^(0), (4) 
where T and V are first order differential operators defined by 

Their exponentiations, e eT and e eV , are displacement operators which displace and p^ 
forward in time via 

q% -> ?i + epi and p» — > p» + eFj. (6) 

Thus, if e^ T+v ^ can be factorized into products of the displacement operators e eT and e eV , 
each such factorization gives rise to an algorithm for evolving the system forward in time. 
For example, the second order factorization 

T 2 (e) = e ^ T e eV eh T = exp[e(T + V) + e 3 C + +0(e 5 ) ■ ■ •], (7) 

corresponds to the second order algorithm 

1 

qi = Qo + 2 e Po> 
pi = p + eF(qi), 

q2 = qi + 2 e Pi> ( 8 ) 

where qo, po and q 2 , pi are the initial and final states of the algorithm respectively. This 
second order symplectic algorithm only requires one evaluation of the force. 

The bilaterally symmetric form of ^(e) automatically guarantees that it is time- 
reversible, 



T 2 (-e)T 2 (e) = 1, (9) 

and implies that log^) can only be an odd function of e, as indicated in ([!]). The explicit 
form of the operator C is not needed for our present discussion. 
Consider now the the symmetric triple product 

T 2 (5)T 2 (-s5)T 2 (5) = exp[(2 - s)S(T + V) + (2 - s 3 )5 3 C + 0(5 5 ) + ■■■]. (10) 

This algorithm evolves the system forward for time 5, backward for time sS and forward again 
for time 5. Since it is manifestly time-reversible, its error terms must be odd powers of 5 only 
Morever, its leading first and third order terms can only be the sum of the first and third 
order terms of each constituent algorithm as indicated. This is because non-additive terms 
must come from commutators of operators and the lowest order non-vanishing commutator 
has to have two first order terms and one third order term, which is fifth order. The form of 
(P~0|) naturally suggests that the third order error term can be made to vanish by choosing 

S = 2 1 / 3 . (11) 

Thus if we now rescale S back to the standard step size by setting e = (2 — s)5, the resulting 
triplet product would be correct to 4th order, 

T 4 = T 2 (-^-)tJ^-)T 2 (-^-) = exp[e(T + V) + 0(e 5 ) + ■■■]. (12) 



>2-s' V2- s / "V2-s 
Expanding out the T 2 's gives the explicit form: 

T 4 = e ^T ebl eV e a 2 eT eb2 eV e a 2 eT ebl eV eai eT (13) 

where, by inspection 

1 1 1S ~ 1 U 1 A U S flA\ 

a± = , a 2 = , b\ = , and o 2 = . (14) 

22-s' 22-s' 2-s 2-s v ' 

This fourth order symplectic algorithm was apparently obtained by E. Forest in 1987. How- 



ever, its original derivation was very complicated and was not published with Ruth JTT 



until 1990. During this period many groups, including Campostrini and Rossi |12| in 1990, 



Candy and Rozmous |T3[ in 1991, independently published the same algorithm. Our dis- 
cussion followed the earliest published derivation of this algorithm by Creutz and Gocksch 
H in 1989. After they were informed of this algorithm by Campostrini, they provided the 
triplet construction and generalized it to higher order. The triplet construction was also 
independently published by Suzuki [7| and Yoshida || in 1990. 

Higher order algorithms can be obtained by repeating this construction. Starting with 
any nth order symmetric algorithm, 

T n (e) = exp[e(T + V) + e^D + •••], (15) 

the triplet product 

T n (5)T n (-s5)T n (5) = exp[(2 - s)5(T + V) + (2- s n+1 )5 n+1 D + 0(5 n+3 ) + •••], (16) 
will be of order (n + 2) if we choose 

s = 2 1 /(-+i) ( 17 ) 
and renormalize 5 = e/(2 — s) as before. 



III. FORCE GRADIENT ALGORITHMS 

The method of operator factorization can be applied to many different classes of evolution 
equations. However, the triplet concatenations with a negative time step are a special 
construction with more limited applicability. For example, one cannot use it to derive similar 
Diffusion Monte Carlo or finite temperature path integral algorithms, because one cannot 
simulate diffusion backward in time nor sample configurations with negative temperatures. 
The triplet construction is a special example of Suzuki's |H| general proof that, beyond 
second order, it is impossible to factorize e e( - T+ir) only into products of e tT 's and e tV 's 
without introducing negative time steps. For symplectic algorithms this means that one 
can never develop a purely positive time step fourth order algorithm by evaluating only 
the force. For many years the Forest-Ruth (FR) algorithm was the only known fourth 
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order symplectic algorithm. Recently, a deeper understanding of the operator factorization 
process has yielded three new symplectic algorithms || all with purely positive time steps. 
These new algorithms circumvented Suzuki's no-go theorem by evaluating the force and 
its gradient. This corresponds to factorizing e t( - T+v > in terms of operators T, V, and the 
commutator [V, [T, V]]. The latter corresponds to 

[V,[T,V}} = 2F J ^^ = V l \F\ 2 ^, (18) 
dqj dpi Opt 

which is the gradient of the squared magnitude of the force. Of the three algorithms derived 
by Chin ||, algorithm C is particularly outstanding and corresponds to the factorization 

e e(T+V) = e elT e e|y e 4T e eIy e 4T e e|y e elT + (1Q) 

where 

V = V+^e 2 [V,[T,V}]. (20) 
The algorithm itself can be read off directly as 



qi= qo + -epo, 

o 

Pi= Po + gCF(cu) 
1 

q2= qi + ^epi 
l 



l 



p 2 = Pl + -e[F(q 2 ) + _ e 2 V|F(q 2 )| 2 j (21) 
1 

qs= q2 + -ep 2 

3 

Ps= P2 + gCF(q 3 ) 
1 

q4= q3 + ^ep3- 



In it was shown that the maximum energy error for this algorithm, when used to solve 
Kepler's problem, is smaller than that of the FR algorithm by a factor of 80. At the moment 
there is no general method for constructing higher order algorithms with only positive time 
steps. It is not even known whether a positive time step 6th order algorithm exists. Thus, 
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beyond 4th order the triplet contraction is still the only systematic way of generating higher 
order algorithms. In this work we show that intrinsic error functions associated with higher 
order algorithms generated from Chin's algorithm C are far smaller than those generated 
from the FR algorithm. 



IV. THE ENERGY AND THE LRL VECTOR 

We gauge the numerical effectiveness of each algorithm by solving the two dimensional 
Kepler problem 

cPq q 

1? = (22) 
with initial conditions qo = (10, 0) and po = (0, 1/10). The resulting highly eccentric (e=0.9) 
orbit provides a non-trivial testing ground for trajectory integration. 

A symmetric nth order symplectic algorithm evolves this system forward in time with 
Hamiltonian 

H{p, q) = H Q (p, q) + e n H n (p, q) + 0(e n+2 ), (23) 

which deviates from the exact Hamiltonian H (p,q) = ^p 2 — l/|q| by an error term 
e n H n (p, q) as indicated. To gauge the intrinsic merit of each algorithm we compare their 
step-size independent error coefficient H n (p,q). This can be extracted numerically as fol- 
lows. Let's start the system with total energy E = H o (p(0), q(0)). Since the Hamiltonian 
(p3|) is conserved by the algorithm, we have 

E = H (p(t), q(t)) + e n H n (p(t), q(f)) + 0(e™+ 2 ) (24) 

Denoting E{t) = H (p(t), q(i)) and H n (t) = H n (p(t), q(t)), we therefore have 

H n (t) = -hm^[E(t)-E }. (25) 

Energy conservation does not directly measure how well the orbit is determined. When 
the time step is not too small, a very noticeable error is that the orbit precesses. One can, 



but it is tedious, directly monitor this orbital precession ||]. It is more expedient to follow 
the rotation of the Laplace-Runge-Lenz (LRL) vector: 



A = p x L — r. 



(26) 



When the orbit is exact the LRL vector is constant, pointing along the semi-major axis of 
the orbit. When the orbit precesses the LRL vector rotates correspondingly. 
For an nth order algorithm, 



dA / ,9A dH„ dA dH n 



(27) 



(28) 
(29) 



dt z ~ jK dq i dpi dpi dqi 
Thus, the rate of change of each component of the LRL vector is of order e n . The components 
themselves, which are time integrals of the above modulo a constant term, must also be of 
order e n . Let the LRL vector initially be of length A Q and lie along the x-axis, then we have 

A x (t)=A + e n A nx (t) + O(e n+2 ), 
A y (t)=e n A ny (t)+0(e n+2 ). 

Since the square of the LRL vector is related to the energy by 

A 2 = 2L 2 E + 1, (30) 

the longitudinal deviation coefficient A nx (t) is related to the energy error coefficient by 

A nx (t) = - E °) = -^r H n(t), (si) 

which gives no new information. The perpendicular deviation coefficient A ny (t) is best 
measured in terms of the rotation angle: 



8(t) = tan" 1 





= e n 




lA x (t)_ 




I A \ 



+ ■ 



(32) 



To compare algorithms we again extract and compare their rotation error coefficient function 
6 n (t) = A ny (t)/A via 



6 n (t) = hm -6(t). 



(33) 



Since this rotational angle is related to some integral of the energy error function, it is a 
better measure of the overall error of the algorithm. 
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V. RESULTS OF COMPARING HIGHER ORDER ALGORITHMS 



By use of the triplet construction, we generated 6th, 8th, 10th, and 12th order algorithms 
from both the Forest-Ruth and Chin's C algorithm. We computed the fractional energy 
deviation, which is just the negative of the energy error coefficient normalized by the initial 



Smaller and smaller time steps e are used until the extracted coefficient function is stablized 
independent of the time step size. This typically occurs in the neighborhood of e = P/5000, 
where P is the period of the orbit. 

Fig.^ compares the (negative) normalized error coefficient functions for the 4th order 
Runge-Kutta, Forest-Ruth and Chin's C algorithms over one period of the orbit. The error 
function for the two symplectic algorithms are substantial only near mid period when the 
particle is at its closest approach to the attractive center. For symplectic algorithms energy 
is conserved over one period, or its non-conservation is periodic. Its average energy error is 
bounded and constant as a function of time. In contrast, the 4th order Runge-Kutta energy 
error function is an irreversible, step-like function over one period. Each successive period 
will increase the error by the same amount resulting in a linearly rising, staircase-like error 
function in time. As noted earlier, the maximum error in Chin's algorithm C is smaller than 
that of the FR algorithm by a factor of 80. However, this error height comparison at one 
point is not meaningful. It is better to compare the energy error averaged over one period. 
This would require the integral of the energy error function. On this basis algorithm C will 
be better still. While the energy error integral can be done, the same goal can be achieved 
by monitoring the rotation of the LRL vector. 

Fig.0 shows the corresponding error coefficient functions of the rotational angle of the 
LRL vector. After each period, the algorithms rotate the LRL vector by a definite amount. 
The error coefficient provides an intrinsic, step-size independent way of comparing this 



energy, 




- 1 



H n (t) 
E 



(34) 
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rotation. In Fig.|2| the rotated angle produced by algorithm C is too small to be visible when 
plotted on the same scale as the other algorithms. The insert gives an enlargement of the 
details. The rotational angle of the LRL vector appears to be related to some integral of 
energy error function. Although we have not been able to demonstrate this analytically, 
numerical integration of the energy error function does give a function similiar in shape 
to the angle coefficient function, having the same numbers of maxima and minima. For 
the Runge-Kutta, Forest-Ruth and Chin's C 4th order algorithms, the magnitudes of this 
rotation cofficient after one period are 2.666, 10.860, and 0.004 respectively. On this basis, 
algorithm C is better than FR by a factor of ~ 3000. When the orbit is integrated over many 
periods the rotational angle from symplectic algorithms increases linearly in a staircase- like 
manner with time. In contrast, the rotational angle of the Runge-Kutta algorithm shows 
a quadratic increase over long times, such as a few thousand periods. This result is easy 
to understand if the rotational angle is related to some integral of the energy error. This 
quadratic increase in the rotation angle of the LRL vector clearly mirrors the quadratic 
increase in phase error of the Runge-Kutta algorithm, as discussed by Gladman, Duncan 
and Candy |§. 

Fig.^ and Fig.[| show the results when both the Forest-Ruth and Chin's C algorithms are 
iterated to 6th order by the triplet product construction. Inserts in both detail algorithm 
C's intricate structure. As an added comparison we also included results for Yoshida's || 
6th order algorithm A, which is a product of 7 second order algorithms (H), some with 
negative time steps. For Yoshida's algorithm, Forest-Ruth and Chin's C algorithm in 6th 
order, the magnitudes of the rotation cofhcients after one period are 11.44, 335.1, and 0.1156 
respectively. Yoshida's algorithm is a factor of 30 better than FR, but algorithm C is a factor 
of 3000 better. Note that if the energy error function is related to the differential of the the 
angle error function, the zeros of the former would correspond to the extrema of the latter. 
The four zero crossings of algorithm C's energy error function are clearly reflected in the 
two maxima and two minima of the corresponding angle error function. 

Fig.[| and Fig.^ give results for the 8th order iterated algorithms based on the Forest- 
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Ruth and Chin's C algorithm. The magnitudes of the angle error coefficients are 1.386 x 10 4 
and 0.4532 respectively, giving a ratio of approximately 3 x 10 4 . Algorithm C retains its 
characteristic shape in both the energy and the angle error function. Fig.[7] and Fig.[8] give 
the corresponding results for the iterated 10th order algorithms. Here the intricate structure 
in the C algorithm is beginning to be washed out. At this high order, quadruple numeric 
precision is necessary to extract these coefficient functions smoothly. The magnitudes of 
the angle error coefficients are now 7.141 x 10 5 and 17.89 respectively, giving a ratio of 
4 x 10 4 . Fig.p and Fig.|l0| give similar results for the 12th order algorithms. At this point all 
structures in the C algorithm are gone. The magnitudes of the angle error coefficients are 
now 4.473 x 10 7 and 427.5 respectively, giving a ratio of 1.05 x 10 5 . 

The iteration of algorithms A and B of Chin || also produced results that are better 
than FR based algorithms. However, we do not detail their results here because they are at 
least one or two orders of magnitude inferior to algorithm C. 



VI. CONCULSIONS 

In this work we have shown that higher order force gradient symplectic algorithms appear 
to be superior to non-gradient symplectic alogorithms as measured by eneregy conservation 
and the rotation of the LRL vector. While it has been shown earlier that 4th order force 
gradient algorithms have smaller energy error coefficients f9|, it was not known that this 
advantage would mulitply dramatically when algorithms are iterated to higher orders. The 
conclusion that one should draw may not be that force gradient algorithms are better, but 
that higher order non-gradient algorithms are far from optimal. Secondly, we suggested that 
the rotation of the LRL vector gives an intergrated measure of an algorithm's merit when 
tested on the Kepler problem. 

The high accuracy of this class of algorithms seemed ideal for long time integration of 
few-body problems, such as that of the solar system ||. For such few-body problems, the 
evaluation of the force gradient is not excessively difficult. It would be useful to examine 
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the merit of this class of algorithms in more physical applications. The distinct advantage 
uncovered in this work, that it is better to iterate a 4th order algorithm with all positive 
time steps, gives further impetus to search for an all positive time step 6th order symplectic 
algorithm. 
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FIG. 1. The normalized energy deviation of a particle in a Keplerian orbit, which measures the 
step-size independent energy error coefficient — i?4(p(i), q(i))/£b- P is the period of the elliptical 
orbit and e is the time step size. RK4, FR, and C denote results for the 4th order Runge-Kutta, 
Forest-Ruth, and Chin's C algorithm respectively. The maximum deviations for algorithm FR and 
C are 21 and 0.27 respectively. 



14 




FIG. 2. The step-size independent error coefficient of the rotation angle of the 
Laplace-Rungc-Lenz vector for 4th order algorithms. The LRL vector rotates substantial only 
when the particle is near mid period, closest to the attractive center. The insert make visible the 
fine structure produced by algorithm C. 
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FIG. 3. The normalized energy deviation for 6th order algorithms. RF and C denote 6th order 
algorithms generated by a triplet product of corresponding 4th algorithms in Fig.Q. The insert 
makes visible the energy devivation of algorithm C, which is not visible in the bigger graph. The 
maximum deviations for algorithms FR, Y and C are 513, 13.6, and 0.74 respectively. 
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FIG. 4. The step-size independent error coefficient of the rotation angle of the LRL vector for 
6th order algorithms as described in Fig.||[ The insert makes visible the minute rotation coefficient 
produced by algorithm C. 
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FIG. 5. The normalized energy deviations for 8th order algorithms, as generated by a triplet 
product of 6th order algorithm described in Fig.||. The insert makes visible the minute energy 
deviation of algorithm C. 
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FIG. 6. The step-size independent error coefficient of the rotation angle of the LRL vector for 
8th order algorithms as described in Fig|| In this order, the algorithm C based algorithm begins 
to rotate in the same sense as the FR base algorithm. 
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FIG. 7. The normalized energy deviations for 10th order algorithms, as generated by a triplet 
product of 8th order algorithm described in FigJ||. The insert shows that the characteristic oscil- 
lations of algorithm C are beginning to disappear. 
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FIG. 8. The step-size independent error coefficient of the rotation angle of the LRL vector for 
10th order algorithms as described in Fig.0. The insert shows that the error coefficient of algorithm 
C begins to look like that of the FR algorithm. 
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FIG. 9. The normalized energy deviations for 12th order algorithms, as generated by a triplet 
product of 10th order algorithm described in Fig.0. The insert shows that there is no longer any 
distinctive structure produced by algorithm C. 
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FIG. 10. The step-size independent error coefficient of the rotation angle of the LRL vector for 
12th order algorithms as described in Fig.||. The insert shows that both algorithms have converged 
to a similar step-like error function. 
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